library(tidyverse)
library(skimr)
library(brms)
style <- read_csv("../input/Styleelongation2017_ZTslh.csv") %>%
complete(day, floret, treatment, time)
── Column specification ───────────────────────────────────────────────────────────────────────────────────
cols(
day = col_double(),
floret = col_double(),
treatment = col_character(),
time = col_time(format = ""),
`length (mm)` = col_double(),
ZT = col_double()
)
style <- style %>% rename(trt=treatment, length=`length (mm)`) %>%
mutate(hour=as.numeric(time-min(time)) / 3600, # elapsed hours, first hour = 0
id=str_c(floret,day,trt,sep="-")) # unique id for each floret
style %>% mutate(trt=as.factor(trt)) %>% skim()
── Data Summary ────────────────────────
Values
Name Piped data
Number of rows 1368
Number of columns 8
_______________________
Column type frequency:
character 1
difftime 1
factor 1
numeric 5
________________________
Group variables None
── Variable type: character ───────────────────────────────────────────────────────────────────────────────
skim_variable n_missing complete_rate min max empty n_unique whitespace
1 id 0 1 8 12 0 72 0
── Variable type: difftime ────────────────────────────────────────────────────────────────────────────────
skim_variable n_missing complete_rate min max median n_unique
1 time 0 1 19800 secs 36000 secs 27900 secs 19
── Variable type: factor ──────────────────────────────────────────────────────────────────────────────────
skim_variable n_missing complete_rate ordered n_unique top_counts
1 trt 0 1 FALSE 3 Eas: 456, Wes: 456, Wes: 456
── Variable type: numeric ─────────────────────────────────────────────────────────────────────────────────
skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
1 day 0 1 5.12 2.71 1 2.75 5.5 7.25 9 ▇▃▃▇▇
2 floret 0 1 2 0.817 1 1 2 3 3 ▇▁▇▁▇
3 length 45 0.967 8.47 1.37 6.24 7.40 8.07 9.35 13.7 ▇▇▃▁▁
4 ZT 39 0.971 1.19 1.34 -1 0 1.25 2.25 3.5 ▇▇▆▇▇
5 hour 0 1 2.25 1.37 0 1 2.25 3.5 4.5 ▇▇▆▇▇
unique(style$trt)
[1] "East" "West" "WestHeat"
Plot to see what we have:
style %>%
ggplot(aes(x=time, y=length, color=trt, shape=as.character(floret))) +
geom_point() +
facet_wrap(~day)
Warning: Removed 45 rows containing missing values (geom_point).
remove a few missing data rows that are not at the end of the series
style <- style %>% filter(day!=1 | floret!=3 | trt!="East" | hour!=3.5) %>%
filter(!(is.na(length) & hour < 2))
Fix height to max height. This is necessary because our growth model can’t account for shrinkage after reaching max height. Fill in missing end-of-series values, but note that they are censored.
style <- style %>%
mutate(censored = ifelse(is.na(length), "right", "none")) %>%
group_by(id) %>%
arrange(hour) %>%
mutate(length2=ifelse((length < max(length, na.rm = TRUE) | is.na(length)) & row_number() > which.max(length),
max(length, na.rm = TRUE), length)) %>%
ungroup() %>%
filter(!is.na(length2))
Plot smoothed data
style %>%
ggplot(aes(x=time, y=length2, color=trt)) +
geom_smooth() +
facet_wrap(~day)
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
And plot averaged over days:
style %>%
ggplot(aes(x=time, y=length2, color=trt)) +
geom_smooth()
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
To help figure out reasonable priors
weibull <- function(Lmin, Lmax, k, delta, hour) {
y <- Lmax - (Lmax - Lmin) * exp(-(k*hour)^delta)
tibble(hour=hour, y=y)
}
set up tibble with some priors
x <- seq(0,10,.01)
plot(x, dexp(x,.1), type = "l")
k only
tibble(
Lmin=7.5,#rnorm(50, 7.5, 1),
Lmax=11.5,#rnorm(50, 11.5, 1),
k=rexp(50,1),
delta=1) %>% #rnorm(50,1,.5)) %>%
mutate(delta=ifelse(delta < 0, 0, delta)) %>%
mutate(prediction=pmap(list(Lmin, Lmax, k, delta), weibull, seq(0,4.5,.1) ),
sample=row_number()) %>%
unnest(prediction) %>%
ggplot(aes(x=hour,y=y, group=sample)) +
geom_line(alpha=.2)
delta only
tibble(
Lmin=7.5,#rnorm(50, 7.5, 1),
Lmax=11.5,#rnorm(50, 11.5, 1),
k=.5,
delta=rexp(50,.5)) %>%
#mutate(delta=ifelse(delta < 0, 0, delta)) %>%
mutate(prediction=pmap(list(Lmin, Lmax, k, delta), weibull, seq(0,4.5,.1) ),
sample=row_number()) %>%
unnest(prediction) %>%
ggplot(aes(x=hour,y=y, group=sample)) +
geom_line(alpha=.2)
delta and k
tibble(
Lmin=7.5,#rnorm(50, 7.5, 1),
Lmax=11.5,#rnorm(50, 11.5, 1),
k=rexp(100, 1),
delta=rexp(100,.5)) %>%
mutate(delta=ifelse(delta < 0, 0, delta)) %>%
mutate(prediction=pmap(list(Lmin, Lmax, k, delta), weibull, seq(0,4.5,.1) ),
sample=row_number()) %>%
unnest(prediction) %>%
ggplot(aes(x=hour,y=y, group=sample)) +
geom_line(alpha=.2)
all together
tibble(
Lmin=rnorm(100, 7.5, .25),
Lmax=rnorm(100, 11.5, .25),
k=rexp(100, 1),
delta=rexp(100, 0.5)) %>%
mutate(prediction=pmap(list(Lmin, Lmax, k, delta), weibull, seq(0,4.5,.1) ),
sample=row_number()) %>%
unnest(prediction) %>%
ggplot(aes(x=hour,y=y, group=sample)) +
geom_line(alpha=.2)
what priors are available to be set?
get_prior(bf(
length2 | cens(censored) ~ Lmax - (Lmax - Lmin) * exp(-(k*hour)^delta), #Weibull model
Lmax + Lmin ~ 1,
delta + k ~ trt + (1|id),
nl=TRUE ),
data=style)
prior class coef group resp dpar nlpar bound source
student_t(3, 0, 2.5) sigma default
(flat) b delta default
(flat) b Intercept delta (vectorized)
(flat) b trtWest delta (vectorized)
(flat) b trtWestHeat delta (vectorized)
student_t(3, 0, 2.5) sd delta default
student_t(3, 0, 2.5) sd id delta (vectorized)
student_t(3, 0, 2.5) sd Intercept id delta (vectorized)
(flat) b k default
(flat) b Intercept k (vectorized)
(flat) b trtWest k (vectorized)
(flat) b trtWestHeat k (vectorized)
student_t(3, 0, 2.5) sd k default
student_t(3, 0, 2.5) sd id k (vectorized)
student_t(3, 0, 2.5) sd Intercept id k (vectorized)
(flat) b Lmax default
(flat) b Intercept Lmax (vectorized)
(flat) b Lmin default
(flat) b Intercept Lmin (vectorized)
prior1 <- c(prior(normal(7.5,.25), nlpar="Lmin"),
prior(normal(11.5, .25),nlpar="Lmax"),
prior(exponential(1),nlpar="k", lb=0),
prior(exponential(.5),nlpar="delta", lb=0) # 1 is exponential model
)
fit1c <- brm(
formula=bf(
length2 | cens(censored) ~ Lmax - (Lmax - Lmin) * exp(-(k*hour)^delta), #Weibull model
Lmax + Lmin + k + delta ~ 1, #minimal model, parameters do not vary for treatment or random effects
nl=TRUE),
data=style,
prior=prior1,
sample_prior = "yes",
cores=4,
chains=4)
Compiling Stan program...
Start sampling
starting worker pid=34485 on localhost:11685 at 15:39:20.689
starting worker pid=34499 on localhost:11685 at 15:39:20.983
starting worker pid=34513 on localhost:11685 at 15:39:21.255
starting worker pid=34527 on localhost:11685 at 15:39:21.518
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
Chain 1: Stan can't start sampling from this initial value.
Chain 1: Rejecting initial value:
Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
Chain 1: Stan can't start sampling from this initial value.
Chain 1: Rejecting initial value:
Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
Chain 1: Stan can't start sampling from this initial value.
Chain 1: Rejecting initial value:
Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
Chain 1: Stan can't start sampling from this initial value.
Chain 1: Rejecting initial value:
Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
Chain 1: Stan can't start sampling from this initial value.
Chain 1: Rejecting initial value:
Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
Chain 1: Stan can't start sampling from this initial value.
Chain 1: Rejecting initial value:
Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
Chain 1: Stan can't start sampling from this initial value.
Chain 1: Rejecting initial value:
Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
Chain 1: Stan can't start sampling from this initial value.
Chain 1:
Chain 1: Gradient evaluation took 0.001376 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 13.76 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 0.00087 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 8.7 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3: Rejecting initial value:
Chain 3: Log probability evaluates to log(0), i.e. negative infinity.
Chain 3: Stan can't start sampling from this initial value.
Chain 3:
Chain 3: Gradient evaluation took 0.000872 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 8.72 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 0.000855 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 8.55 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 27.675 seconds (Warm-up)
Chain 2: 28.128 seconds (Sampling)
Chain 2: 55.803 seconds (Total)
Chain 2:
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 29.426 seconds (Warm-up)
Chain 1: 29.167 seconds (Sampling)
Chain 1: 58.593 seconds (Total)
Chain 1:
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 28.75 seconds (Warm-up)
Chain 3: 28.803 seconds (Sampling)
Chain 3: 57.553 seconds (Total)
Chain 3:
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 27.737 seconds (Warm-up)
Chain 4: 30.294 seconds (Sampling)
Chain 4: 58.031 seconds (Total)
Chain 4:
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 5383764 287.6 8292668 442.9 NA 8292668 442.9
Vcells 16841369 128.5 83330576 635.8 16384 203442360 1552.2
fit1c <- add_criterion(fit1c, "loo")
summary(fit1c)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: length2 | cens(censored) ~ Lmax - (Lmax - Lmin) * exp(-(k * hour)^delta)
Lmax ~ 1
Lmin ~ 1
k ~ 1
delta ~ 1
Data: style (Number of observations: 1362)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Lmax_Intercept 11.34 0.23 10.91 11.80 1.00 1487 2019
Lmin_Intercept 7.40 0.05 7.29 7.51 1.00 2226 2007
k_Intercept 0.28 0.01 0.26 0.30 1.00 1596 2264
delta_Intercept 3.00 0.24 2.59 3.53 1.00 1646 1997
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 1.04 0.02 1.01 1.09 1.00 2793 2325
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 5383518 287.6 8292668 442.9 NA 8292668 442.9
Vcells 16603765 126.7 66664461 508.7 16384 203442360 1552.2
plot(fit1c)
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 5501005 293.8 8292668 442.9 NA 8292668 442.9
Vcells 17411913 132.9 53331569 406.9 16384 203442360 1552.2
pairs(fit1c)
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 5676939 303.2 8292668 442.9 NA 8292668 442.9
Vcells 18319951 139.8 53331569 406.9 16384 203442360 1552.2
pred1c <- predict(fit1c)
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 5388768 287.8 8292668 442.9 NA 8292668 442.9
Vcells 16597600 126.7 61578767 469.9 16384 203442360 1552.2
pred1c <- cbind(pred1c, style)
pred1c %>% ggplot(aes(x=hour, y=length2, color=trt)) +
geom_smooth(span=1) +
geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_line(aes(y=Estimate),color="yellow")
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
pp_check(fit1c)
Using 10 posterior samples for ppc type 'dens_overlay' by default.
Warning: Censored responses are not shown in 'pp_check'.
prior2 <- c(prior(normal(7.5,.25), nlpar="Lmin"),
prior(normal(11.5, .25),nlpar="Lmax"),
prior(exponential(1),nlpar="k", lb=0),
prior(exponential(0.5),nlpar="delta"), # 1 is exponential model
prior(normal(0,0.5), nlpar="delta", coef="trtWest"),
prior(normal(0,0.5), nlpar="delta", coef="trtWestHeat")
)
fit2c <- brm( # ignore the warning about lower bounds. If we set it for the intercept, it impacts the other delta coefficents as well and we don't want that.
formula=bf(
length2 | cens(censored) ~ Lmax - (Lmax - Lmin) * exp(-(k*hour)^delta), #Weibull model
Lmax + Lmin + k ~ 1,
delta ~ trt + (1|id), #
nl=TRUE ),
data=style,
prior=prior2,
sample_prior = "yes",
cores=4,
chains=4,
iter=5000 #,
#control=list(adapt_delta=0.9)
)
Warning: It appears as if you have specified a lower bounded prior on a parameter that has no natural lower bound.
If this is really what you want, please specify argument 'lb' of 'set_prior' appropriately.
Warning occurred for prior
b_delta ~ exponential(0.5)
Compiling Stan program...
Start sampling
starting worker pid=35268 on localhost:11685 at 16:11:17.561
starting worker pid=35282 on localhost:11685 at 16:11:17.823
starting worker pid=35296 on localhost:11685 at 16:11:18.135
starting worker pid=35310 on localhost:11685 at 16:11:18.474
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: exponential_lpdf: Random variable is -1.33308, but must be nonnegative! (in 'anon_model', line 78, column 2 to column 47)
Chain 1: Rejecting initial value:
Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
Chain 1: Stan can't start sampling from this initial value.
Chain 1:
Chain 1: Gradient evaluation took 0.002339 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 23.39 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: exponential_lpdf: Random variable is -0.020832, but must be nonnegative! (in 'anon_model', line 78, column 2 to column 47)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: exponential_lpdf: Random variable is -0.763222, but must be nonnegative! (in 'anon_model', line 78, column 2 to column 47)
Chain 2: Rejecting initial value:
Chain 2: Log probability evaluates to log(0), i.e. negative infinity.
Chain 2: Stan can't start sampling from this initial value.
Chain 2: Rejecting initial value:
Chain 2: Log probability evaluates to log(0), i.e. negative infinity.
Chain 2: Stan can't start sampling from this initial value.
Chain 2: Rejecting initial value:
Chain 2: Log probability evaluates to log(0), i.e. negative infinity.
Chain 2: Stan can't start sampling from this initial value.
Chain 2: Rejecting initial value:
Chain 2: Log probability evaluates to log(0), i.e. negative infinity.
Chain 2: Stan can't start sampling from this initial value.
Chain 2: Rejecting initial value:
Chain 2: Log probability evaluates to log(0), i.e. negative infinity.
Chain 2: Stan can't start sampling from this initial value.
Chain 2:
Chain 2: Gradient evaluation took 0.000967 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 9.67 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: exponential_lpdf: Random variable is -1.0895, but must be nonnegative! (in 'anon_model', line 78, column 2 to column 47)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: exponential_lpdf: Random variable is -1.71041, but must be nonnegative! (in 'anon_model', line 78, column 2 to column 47)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: exponential_lpdf: Random variable is -0.560105, but must be nonnegative! (in 'anon_model', line 78, column 2 to column 47)
Chain 3: Rejecting initial value:
Chain 3: Log probability evaluates to log(0), i.e. negative infinity.
Chain 3: Stan can't start sampling from this initial value.
Chain 3:
Chain 3: Gradient evaluation took 0.000778 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 7.78 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: exponential_lpdf: Random variable is -1.3763, but must be nonnegative! (in 'anon_model', line 78, column 2 to column 47)
Chain 4: Rejecting initial value:
Chain 4: Log probability evaluates to log(0), i.e. negative infinity.
Chain 4: Stan can't start sampling from this initial value.
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: exponential_lpdf: Random variable is -1.98786, but must be nonnegative! (in 'anon_model', line 78, column 2 to column 47)
Chain 4: Rejecting initial value:
Chain 4: Log probability evaluates to log(0), i.e. negative infinity.
Chain 4: Stan can't start sampling from this initial value.
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: exponential_lpdf: Random variable is -0.120339, but must be nonnegative! (in 'anon_model', line 78, column 2 to column 47)
Chain 4: Rejecting initial value:
Chain 4: Log probability evaluates to log(0), i.e. negative infinity.
Chain 4: Stan can't start sampling from this initial value.
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: exponential_lpdf: Random variable is -0.883203, but must be nonnegative! (in 'anon_model', line 78, column 2 to column 47)
Chain 4: Rejecting initial value:
Chain 4: Log probability evaluates to log(0), i.e. negative infinity.
Chain 4: Stan can't start sampling from this initial value.
Chain 4:
Chain 4: Gradient evaluation took 0.000991 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 9.91 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 5000 [ 0%] (Warmup)
Chain 1: Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 3: Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 4: Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 1: Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 3: Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 4: Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 3: Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 1: Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 4: Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 2: Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 3: Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 1: Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 4: Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 3: Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 3: Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 1: Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 1: Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 2: Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 4: Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 4: Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 3: Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 1: Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 2: Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 4: Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 3: Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 1: Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 2: Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 2: Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 4: Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 3: Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 1: Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 2: Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 4: Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 3: Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 1: Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 2: Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 4: Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 3: Iteration: 5000 / 5000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 239.44 seconds (Warm-up)
Chain 3: 185.642 seconds (Sampling)
Chain 3: 425.082 seconds (Total)
Chain 3:
Chain 1: Iteration: 5000 / 5000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 248.287 seconds (Warm-up)
Chain 1: 182.865 seconds (Sampling)
Chain 1: 431.152 seconds (Total)
Chain 1:
Chain 2: Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 4: Iteration: 5000 / 5000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 274.278 seconds (Warm-up)
Chain 4: 176.014 seconds (Sampling)
Chain 4: 450.292 seconds (Total)
Chain 4:
Chain 2: Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 2: Iteration: 5000 / 5000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 328.205 seconds (Warm-up)
Chain 2: 146.374 seconds (Sampling)
Chain 2: 474.579 seconds (Total)
Chain 2:
fit2c <- add_criterion(fit2c, "loo")
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 5603569 299.3 8292668 442.9 NA 8292668 442.9
Vcells 17212243 131.4 124350110 948.8 16384 203442360 1552.2
summary(fit2c)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: length2 | cens(censored) ~ Lmax - (Lmax - Lmin) * exp(-(k * hour)^delta)
Lmax ~ 1
Lmin ~ 1
k ~ 1
delta ~ trt + (1 | id)
Data: style (Number of observations: 1362)
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
total post-warmup samples = 10000
Group-Level Effects:
~id (Number of levels: 72)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(delta_Intercept) 1.40 0.16 1.13 1.75 1.00 1074 2637
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Lmax_Intercept 15.89 0.16 15.59 16.20 1.00 6236 6449
Lmin_Intercept 7.31 0.03 7.26 7.36 1.00 9392 7973
k_Intercept 0.18 0.00 0.18 0.18 1.00 5063 6317
delta_Intercept 2.62 0.26 2.14 3.14 1.01 539 1536
delta_trtWest 1.16 0.31 0.56 1.76 1.00 850 2356
delta_trtWestHeat -0.03 0.31 -0.65 0.57 1.02 677 1259
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.55 0.01 0.53 0.58 1.00 8127 7283
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 5603631 299.3 8292668 442.9 NA 8292668 442.9
Vcells 17212454 131.4 99480088 759.0 16384 203442360 1552.2
plot(fit2c, ask = FALSE)
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 5338178 285.1 8292668 442.9 NA 8292668 442.9
Vcells 16153074 123.3 72548821 553.6 16384 171505768 1308.5
plot(conditional_effects(fit2c), ask=FALSE)
pp_check(fit2c)
Using 10 posterior samples for ppc type 'dens_overlay' by default.
Warning: Censored responses are not shown in 'pp_check'.
pred2c <- predict(fit2c)
pred2c %>%
cbind(style) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_smooth(span=1) +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_smooth(aes(y=Estimate), lty=2, span=1)
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
pred2c %>%
cbind(style) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_smooth(span=1) +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_smooth(aes(y=Estimate), lty=2, span=1) +
facet_wrap(~day, scales = "free_y")
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
pred2c %>%
cbind(style) %>%
group_by(day, trt, hour) %>%
summarize(length2 = mean(length2),
Estimate = mean(Estimate)) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_point() +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_line(aes(y=Estimate), lty=2) +
facet_wrap(~day, scales = "free_y")
`summarise()` has grouped output by 'day', 'trt'. You can override using the `.groups` argument.
loo_compare(fit1c, fit2c)
elpd_diff se_diff
fit2c 0.0 0.0
fit1c -812.7 27.7
Fit2c is preferred, but predictions are not that good relative to observed
hypothesis(fit2c, c("delta_trtWest = 0",
"delta_trtWestHeat = 0",
"delta_trtWest = delta_trtWestHeat"))
Hypothesis Tests for class b:
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.
#inits3 <- list(b_Lmax=11.5, b_Lmin=7.5, b_delta=2, b_k=c(.1,0,0), z_1=array(0, dim=c(72,1362)))
inits3 <- function() { list(
b_Lmax=rnorm(1, 11.5, .25),
b_Lmin=rnorm(1, 7.5, .25),
b_delta=rexp(1, .5),
b_k=c(rexp(1, 1), 0, 0),
z_1=array(0, dim=c(72,1362))
)
}
#inits3 <- replicate(4, inits3, simplify = FALSE)
prior3 <- c(prior(normal(7.5,.25), nlpar="Lmin"),
prior(normal(11.5, .25),nlpar="Lmax"),
prior(exponential(1),nlpar="k"),
prior(normal(0,0.1), nlpar="k", coef="trtWest"),
prior(normal(0,0.1), nlpar="k", coef="trtWestHeat"),
prior(exponential(0.5),nlpar="delta")
)
fit3c <- brm( # ignore the warning about k.
formula=bf(
length2 | cens(censored) ~ Lmax - (Lmax - Lmin) * exp(-(k*hour)^delta), #Weibull model
Lmax + Lmin + delta ~ 1,
k ~ trt + (1|id), #
nl=TRUE ),
data=style,
prior=prior3,
sample_prior = "yes",
inits = inits3,
cores=4,
chains=4,
iter=5000,
control=list(adapt_delta=0.99)
)
Warning: It appears as if you have specified a lower bounded prior on a parameter that has no natural lower bound.
If this is really what you want, please specify argument 'lb' of 'set_prior' appropriately.
Warning occurred for prior
b_k ~ exponential(1)
b_delta ~ exponential(0.5)
Compiling Stan program...
Start sampling
starting worker pid=39124 on localhost:11685 at 18:23:17.017
starting worker pid=39138 on localhost:11685 at 18:23:17.253
starting worker pid=39152 on localhost:11685 at 18:23:17.492
starting worker pid=39166 on localhost:11685 at 18:23:17.727
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.000706 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 7.06 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 0.012113 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 121.13 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.000771 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 7.71 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 0.000995 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 9.95 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 5000 [ 0%] (Warmup)
Chain 3: Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 2: Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 4: Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 3: Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 1: Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 2: Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 4: Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 2: Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 3: Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 1: Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 4: Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 2: Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 1: Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 3: Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 4: Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 2: Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 2: Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 1: Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 4: Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 4: Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 3: Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 3: Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 2: Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 3: Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 4: Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 1: Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 1: Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 3: Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 2: Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 4: Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 1: Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 3: Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 3: Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 1: Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 4: Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 2: Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 3: Iteration: 5000 / 5000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 391.743 seconds (Warm-up)
Chain 3: 195.878 seconds (Sampling)
Chain 3: 587.621 seconds (Total)
Chain 3:
Chain 1: Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 4: Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 2: Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 1: Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 4: Iteration: 5000 / 5000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 386.775 seconds (Warm-up)
Chain 4: 259.497 seconds (Sampling)
Chain 4: 646.272 seconds (Total)
Chain 4:
Chain 1: Iteration: 5000 / 5000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 467.077 seconds (Warm-up)
Chain 1: 203.198 seconds (Sampling)
Chain 1: 670.275 seconds (Total)
Chain 1:
Chain 2: Iteration: 5000 / 5000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 363.982 seconds (Warm-up)
Chain 2: 305.675 seconds (Sampling)
Chain 2: 669.657 seconds (Total)
Chain 2:
Warning: There were 13 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess
fit3c <- add_criterion(fit3c, "loo")
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 5801529 309.9 9991201 533.6 NA 9991201 533.6
Vcells 19434076 148.3 125390965 956.7 16384 203442360 1552.2
summary(fit3c)
Warning: There were 13 divergent transitions after warmup. Increasing adapt_delta above 0.99 may help. See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
Family: gaussian
Links: mu = identity; sigma = identity
Formula: length2 | cens(censored) ~ Lmax - (Lmax - Lmin) * exp(-(k * hour)^delta)
Lmax ~ 1
Lmin ~ 1
delta ~ 1
k ~ trt + (1 | id)
Data: style (Number of observations: 1362)
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
total post-warmup samples = 10000
Group-Level Effects:
~id (Number of levels: 72)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(k_Intercept) 0.07 0.01 0.06 0.08 1.01 458 1114
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Lmax_Intercept 13.23 0.12 13.01 13.47 1.00 3324 4655
Lmin_Intercept 7.29 0.03 7.24 7.34 1.00 2989 4500
delta_Intercept 2.51 0.07 2.37 2.65 1.00 2609 4038
k_Intercept 0.26 0.01 0.23 0.29 1.01 270 471
k_trtWest -0.08 0.02 -0.12 -0.04 1.01 298 591
k_trtWestHeat -0.03 0.02 -0.07 0.01 1.00 262 583
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.50 0.01 0.48 0.52 1.00 3692 5473
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 5801593 309.9 9991201 533.6 NA 9991201 533.6
Vcells 19434303 148.3 100312772 765.4 16384 203442360 1552.2
pairs(fit3c, pars="^b_")
plot(fit3c, ask = FALSE)
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 5989463 319.9 9991201 533.6 NA 9991201 533.6
Vcells 21923939 167.3 64200175 489.9 16384 203442360 1552.2
plot(conditional_effects(fit3c), ask=FALSE)
pp_check(fit3c)
Using 10 posterior samples for ppc type 'dens_overlay' by default.
Warning: Censored responses are not shown in 'pp_check'.
pred3c <- predict(fit3c)
pred3c %>%
cbind(style) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_smooth(span=1) +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_smooth(aes(y=Estimate), lty=2, span=1)
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
pred3c %>%
cbind(style) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_smooth(span=1) +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_smooth(aes(y=Estimate), lty=2, span=1) +
facet_wrap(~day, scales = "free_y")
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
pred3c %>%
cbind(style) %>%
group_by(day, trt, hour) %>%
summarize(length2 = mean(length2),
Estimate = mean(Estimate)) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_point() +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_line(aes(y=Estimate), lty=2) +
facet_wrap(~day, scales = "free_y")
`summarise()` has grouped output by 'day', 'trt'. You can override using the `.groups` argument.
loo_compare(fit1c, fit2c, fit3c)
elpd_diff se_diff
fit3c 0.0 0.0
fit2c -153.9 28.1
fit1c -966.7 41.4
inits4 <- function() { list(
b_Lmax=rnorm(1, 11.5, .25),
b_Lmin=rnorm(1, 7.5, .25),
b_delta=c(rexp(1, .5), 0, 0),
b_k=c(rexp(1, 1), 0, 0),
z_1=array(0, dim=c(72,1)),
z_2=array(0, dim=c(72,1))
)
}
prior4 <- c(prior(normal(7.5,.25), nlpar="Lmin"),
prior(normal(11.5, .25),nlpar="Lmax"),
prior(exponential(1),nlpar="k"),
prior(normal(0,0.02), nlpar="k", coef="trtWest"),
prior(normal(0,0.02), nlpar="k", coef="trtWestHeat"),
prior(exponential(0.5), nlpar="delta"), # 1 is exponential model
prior(normal(0,0.25), nlpar="delta", coef="trtWest"),
prior(normal(0,0.25), nlpar="delta", coef="trtWestHeat")
)
fit_4c <- brm(
formula=bf(
length2 | cens(censored) ~ Lmax - (Lmax - Lmin) * exp(-(k*hour)^delta), #Weibull model
Lmax + Lmin ~ 1,
delta + k ~ trt + (1|id), #
nl=TRUE ),
data=style,
prior=prior4,
inits=inits4,
cores=4,
chains=4,
iter=10000, control = list(adapt_delta=0.95)
)
Warning: It appears as if you have specified a lower bounded prior on a parameter that has no natural lower bound.
If this is really what you want, please specify argument 'lb' of 'set_prior' appropriately.
Warning occurred for prior
b_k ~ exponential(1)
b_delta ~ exponential(0.5)
Compiling Stan program...
Start sampling
starting worker pid=43931 on localhost:11685 at 22:42:07.548
starting worker pid=43945 on localhost:11685 at 22:42:07.787
starting worker pid=43959 on localhost:11685 at 22:42:08.018
starting worker pid=43973 on localhost:11685 at 22:42:08.250
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
Chain 1: Stan can't start sampling from this initial value.
Chain 1:
Chain 1: Gradient evaluation took 0.001149 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 11.49 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 10000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 0.001544 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 15.44 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 10000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.001653 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 16.53 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 10000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4: Rejecting initial value:
Chain 4: Log probability evaluates to log(0), i.e. negative infinity.
Chain 4: Stan can't start sampling from this initial value.
Chain 4:
Chain 4: Gradient evaluation took 0.001191 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 11.91 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 10000 [ 0%] (Warmup)
Chain 4: Iteration: 1000 / 10000 [ 10%] (Warmup)
Chain 3: Iteration: 1000 / 10000 [ 10%] (Warmup)
Chain 4: Iteration: 2000 / 10000 [ 20%] (Warmup)
Chain 3: Iteration: 2000 / 10000 [ 20%] (Warmup)
Chain 2: Iteration: 1000 / 10000 [ 10%] (Warmup)
Chain 4: Iteration: 3000 / 10000 [ 30%] (Warmup)
Chain 2: Iteration: 2000 / 10000 [ 20%] (Warmup)
Chain 3: Iteration: 3000 / 10000 [ 30%] (Warmup)
Chain 4: Iteration: 4000 / 10000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 10000 [ 10%] (Warmup)
Chain 2: Iteration: 3000 / 10000 [ 30%] (Warmup)
Chain 3: Iteration: 4000 / 10000 [ 40%] (Warmup)
Chain 4: Iteration: 5000 / 10000 [ 50%] (Warmup)
Chain 4: Iteration: 5001 / 10000 [ 50%] (Sampling)
Chain 2: Iteration: 4000 / 10000 [ 40%] (Warmup)
Chain 3: Iteration: 5000 / 10000 [ 50%] (Warmup)
Chain 3: Iteration: 5001 / 10000 [ 50%] (Sampling)
Chain 4: Iteration: 6000 / 10000 [ 60%] (Sampling)
Chain 2: Iteration: 5000 / 10000 [ 50%] (Warmup)
Chain 2: Iteration: 5001 / 10000 [ 50%] (Sampling)
Chain 3: Iteration: 6000 / 10000 [ 60%] (Sampling)
Chain 4: Iteration: 7000 / 10000 [ 70%] (Sampling)
Chain 2: Iteration: 6000 / 10000 [ 60%] (Sampling)
Chain 3: Iteration: 7000 / 10000 [ 70%] (Sampling)
Chain 4: Iteration: 8000 / 10000 [ 80%] (Sampling)
Chain 2: Iteration: 7000 / 10000 [ 70%] (Sampling)
Chain 3: Iteration: 8000 / 10000 [ 80%] (Sampling)
Chain 4: Iteration: 9000 / 10000 [ 90%] (Sampling)
Chain 2: Iteration: 8000 / 10000 [ 80%] (Sampling)
Chain 3: Iteration: 9000 / 10000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 10000 [ 20%] (Warmup)
Chain 4: Iteration: 10000 / 10000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 718.1 seconds (Warm-up)
Chain 4: 849.172 seconds (Sampling)
Chain 4: 1567.27 seconds (Total)
Chain 4:
Chain 2: Iteration: 9000 / 10000 [ 90%] (Sampling)
Chain 3: Iteration: 10000 / 10000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 817.137 seconds (Warm-up)
Chain 3: 798.559 seconds (Sampling)
Chain 3: 1615.7 seconds (Total)
Chain 3:
Chain 1: Iteration: 3000 / 10000 [ 30%] (Warmup)
Chain 2: Iteration: 10000 / 10000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 963.456 seconds (Warm-up)
Chain 2: 726.01 seconds (Sampling)
Chain 2: 1689.47 seconds (Total)
Chain 2:
Chain 1: Iteration: 4000 / 10000 [ 40%] (Warmup)
Chain 1: Iteration: 5000 / 10000 [ 50%] (Warmup)
Chain 1: Iteration: 5001 / 10000 [ 50%] (Sampling)
Chain 1: Iteration: 6000 / 10000 [ 60%] (Sampling)
Chain 1: Iteration: 7000 / 10000 [ 70%] (Sampling)
Chain 1: Iteration: 8000 / 10000 [ 80%] (Sampling)
Chain 1: Iteration: 9000 / 10000 [ 90%] (Sampling)
Chain 1: Iteration: 10000 / 10000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 1751.66 seconds (Warm-up)
Chain 1: 275.371 seconds (Sampling)
Chain 1: 2027.03 seconds (Total)
Chain 1:
Warning: There were 13 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
fit_4c <- add_criterion(fit_4c, "loo")
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 6038215 322.5 9991201 533.6 NA 9991201 533.6
Vcells 27018499 206.2 239593280 1828.0 16384 352107441 2686.4
summary(fit_4c)
Warning: There were 13 divergent transitions after warmup. Increasing adapt_delta above 0.95 may help. See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
Family: gaussian
Links: mu = identity; sigma = identity
Formula: length2 | cens(censored) ~ Lmax - (Lmax - Lmin) * exp(-(k * hour)^delta)
Lmax ~ 1
Lmin ~ 1
delta ~ trt + (1 | id)
k ~ trt + (1 | id)
Data: style (Number of observations: 1362)
Samples: 4 chains, each with iter = 10000; warmup = 5000; thin = 1;
total post-warmup samples = 20000
Group-Level Effects:
~id (Number of levels: 72)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(delta_Intercept) 0.44 0.08 0.30 0.60 1.00 10160 14357
sd(k_Intercept) 0.07 0.01 0.06 0.08 1.00 6734 10868
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Lmax_Intercept 13.37 0.12 13.14 13.61 1.00 29252 15954
Lmin_Intercept 7.30 0.03 7.25 7.35 1.00 33915 15137
delta_Intercept 2.53 0.12 2.31 2.77 1.00 20251 15535
delta_trtWest 0.22 0.15 -0.07 0.53 1.00 26303 16784
delta_trtWestHeat 0.07 0.14 -0.20 0.34 1.00 22649 16015
k_Intercept 0.24 0.01 0.21 0.26 1.00 4716 8526
k_trtWest -0.04 0.01 -0.06 -0.01 1.00 6706 10191
k_trtWestHeat -0.00 0.01 -0.03 0.02 1.00 6205 9494
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.48 0.01 0.47 0.50 1.00 25540 15145
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 6038280 322.5 9991201 533.6 NA 9991201 533.6
Vcells 27018753 206.2 191674624 1462.4 16384 352107441 2686.4
pairs(fit_4c, pars = "^b_")
plot(fit_4c, ask = FALSE)
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 6193507 330.8 9991201 533.6 NA 9991201 533.6
Vcells 31694194 241.9 153339700 1169.9 16384 352107441 2686.4
plot(conditional_effects(fit_4c), ask=FALSE)
pp_check(fit_4c)
Using 10 posterior samples for ppc type 'dens_overlay' by default.
Warning: Censored responses are not shown in 'pp_check'.
pred_4c <- predict(fit_4c)
pred_4c %>%
cbind(style) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_smooth(span=1) +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_smooth(aes(y=Estimate), lty=2, span=1)
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
pred_4c %>%
cbind(style) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_smooth(span=1) +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_smooth(aes(y=Estimate), lty=2, span=1) +
facet_wrap(~day, scales="free_y")
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
pred_4c %>%
cbind(style) %>%
group_by(trt, hour) %>%
summarize(length2 = mean(length2),
Estimate = mean(Estimate)) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_line() +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_line(aes(y=Estimate), lty=2)
`summarise()` has grouped output by 'trt'. You can override using the `.groups` argument.
pred_4c %>%
cbind(style) %>%
group_by(day, trt, hour) %>%
summarize(length2 = mean(length2),
Estimate = mean(Estimate)) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_point() +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_line(aes(y=Estimate), lty=2) +
facet_wrap(~day, scales = "free_y")
`summarise()` has grouped output by 'day', 'trt'. You can override using the `.groups` argument.
loo_compare(fit1c, fit2c, fit3c, fit_4c)
elpd_diff se_diff
fit_4c 0.0 0.0
fit3c -22.8 6.2
fit2c -176.7 27.8
fit1c -989.5 42.2
inits5 <- function() { list(
b_Lmax=rnorm(1, 11.5, .25),
b_Lmin=rnorm(1, 7.5, .25),
b_delta=c(rexp(1, .5), 0, 0),
b_k=c(rexp(1, 1), 0, 0),
z_1=array(0, dim=c(8,1)),
z_2=array(0, dim=c(72,1)),
z_3=array(0, dim=c(72,1))
)
}
prior5 <- c(prior(normal(7.5,.25), nlpar="Lmin"),
prior(normal(11.5, .25),nlpar="Lmax"),
prior(exponential(1),nlpar="k"),
prior(normal(0,0.02), nlpar="k", coef="trtWest"),
prior(normal(0,0.02), nlpar="k", coef="trtWestHeat"),
prior(exponential(0.5), nlpar="delta"), # 1 is exponential model
prior(normal(0,0.25), nlpar="delta", coef="trtWest"),
prior(normal(0,0.25), nlpar="delta", coef="trtWestHeat")
)
fit_5c <- brm(
formula=bf(
length2 | cens(censored) ~ Lmax - (Lmax - Lmin) * exp(-(k*hour)^delta), #Weibull model
Lmin ~ 1,
Lmax ~ (1|day),
delta + k ~ trt + (1|id), #
nl=TRUE ),
data=style,
prior=prior5,
sample_prior = "yes",
inits=inits5,
cores=4,
chains=4,
iter=2000
)
Warning: It appears as if you have specified a lower bounded prior on a parameter that has no natural lower bound.
If this is really what you want, please specify argument 'lb' of 'set_prior' appropriately.
Warning occurred for prior
b_k ~ exponential(1)
b_delta ~ exponential(0.5)
Compiling Stan program...
Start sampling
starting worker pid=72742 on localhost:11685 at 10:44:09.560
starting worker pid=72756 on localhost:11685 at 10:44:09.797
starting worker pid=72770 on localhost:11685 at 10:44:10.018
starting worker pid=72784 on localhost:11685 at 10:44:10.242
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.001723 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 17.23 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2: Rejecting initial value:
Chain 2: Log probability evaluates to log(0), i.e. negative infinity.
Chain 2: Stan can't start sampling from this initial value.
Chain 2:
Chain 2: Gradient evaluation took 0.001255 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 12.55 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3: Rejecting initial value:
Chain 3: Log probability evaluates to log(0), i.e. negative infinity.
Chain 3: Stan can't start sampling from this initial value.
Chain 3:
Chain 3: Gradient evaluation took 0.003629 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 36.29 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 0.003214 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 32.14 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 223.221 seconds (Warm-up)
Chain 1: 91.087 seconds (Sampling)
Chain 1: 314.308 seconds (Total)
Chain 1:
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 224.296 seconds (Warm-up)
Chain 3: 94.823 seconds (Sampling)
Chain 3: 319.119 seconds (Total)
Chain 3:
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 218.1 seconds (Warm-up)
Chain 4: 98.073 seconds (Sampling)
Chain 4: 316.173 seconds (Total)
Chain 4:
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 678.57 seconds (Warm-up)
Chain 2: 145.349 seconds (Sampling)
Chain 2: 823.919 seconds (Total)
Chain 2:
fit_5c <- add_criterion(fit_5c, "loo")
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 6507990 347.6 9991201 533.6 NA 9991201 533.6
Vcells 32971270 251.6 130397925 994.9 16384 428446769 3268.8
summary(fit_5c)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: length2 | cens(censored) ~ Lmax - (Lmax - Lmin) * exp(-(k * hour)^delta)
Lmin ~ 1
Lmax ~ (1 | day)
delta ~ trt + (1 | id)
k ~ trt + (1 | id)
Data: style (Number of observations: 1362)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Group-Level Effects:
~day (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Lmax_Intercept) 2.03 0.56 1.24 3.39 1.01 604 820
~id (Number of levels: 72)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(delta_Intercept) 1.85 0.28 1.37 2.49 1.00 951 1763
sd(k_Intercept) 0.03 0.00 0.02 0.03 1.00 1032 1923
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Lmin_Intercept 7.36 0.02 7.32 7.40 1.00 2565 2544
Lmax_Intercept 11.47 0.23 11.02 11.94 1.00 1305 2127
delta_Intercept 4.35 0.36 3.70 5.10 1.00 885 1889
delta_trtWest -0.00 0.23 -0.45 0.46 1.00 1892 2752
delta_trtWestHeat 0.06 0.22 -0.39 0.50 1.00 2013 2535
k_Intercept 0.32 0.01 0.31 0.34 1.00 998 1918
k_trtWest -0.09 0.01 -0.10 -0.07 1.00 868 1515
k_trtWestHeat -0.02 0.01 -0.04 -0.01 1.00 922 1379
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.45 0.01 0.43 0.47 1.00 3597 2868
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 6507679 347.6 9991201 533.6 NA 9991201 533.6
Vcells 32666777 249.3 104318340 795.9 16384 428446769 3268.8
plot(fit_5c, ask = FALSE)
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 6693858 357.5 9991201 533.6 NA 9991201 533.6
Vcells 33875577 258.5 104318340 795.9 16384 428446769 3268.8
plot(conditional_effects(fit_5c), ask=FALSE)
pp_check(fit_5c)
Using 10 posterior samples for ppc type 'dens_overlay' by default.
Warning: Censored responses are not shown in 'pp_check'.
pred_5c <- predict(fit_5c)
pred_5c %>%
cbind(style) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_smooth(span=1) +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_smooth(aes(y=Estimate), lty=2, span=1)
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
pred_5c %>%
cbind(style) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_smooth(span=1) +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_smooth(aes(y=Estimate), lty=2, span=1) +
facet_wrap(~day)
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
pred_5c %>%
cbind(style) %>%
group_by(trt, hour) %>%
summarize(length2 = mean(length2),
Estimate = mean(Estimate)) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_line() +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_line(aes(y=Estimate), lty=2)
`summarise()` has grouped output by 'trt'. You can override using the `.groups` argument.
pred_5c %>%
cbind(style) %>%
group_by(day, trt, hour) %>%
summarize(length2 = mean(length2),
Estimate = mean(Estimate)) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_point() +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_line(aes(y=Estimate), lty=2) +
facet_wrap(~day, scales = "free_y")
`summarise()` has grouped output by 'day', 'trt'. You can override using the `.groups` argument.
loo_compare(fit1c, fit2c, fit3c, fit_4c, fit_5c)
elpd_diff se_diff
fit_5c 0.0 0.0
fit_4c -108.8 10.8
fit3c -131.5 13.1
fit2c -285.4 30.6
fit1c -1098.2 43.9
inits6 <- function() { list(
b_Lmax=rnorm(1, 11.5, .25),
b_Lmin=rnorm(1, 7.5, .25),
b_delta=c(rexp(1, .5), 0, 0),
b_k=c(rexp(1, 1), 0, 0),
z_1=array(0, dim=c(8,1)),
z_2=array(0, dim=c(8,1)),
z_3=array(0, dim=c(72,1)),
z_4=array(0, dim=c(72,1))
)
}
prior6 <- c(prior(normal(7.5,.25), nlpar="Lmin"),
prior(normal(11.5, .25),nlpar="Lmax"),
prior(exponential(1),nlpar="k"),
prior(normal(0,0.04), nlpar="k", coef="trtWest"),
prior(normal(0,0.04), nlpar="k", coef="trtWestHeat"),
prior(exponential(0.5), nlpar="delta"), # 1 is exponential model
prior(normal(0,0.25), nlpar="delta", coef="trtWest"),
prior(normal(0,0.25), nlpar="delta", coef="trtWestHeat")
)
fit_6c <- brm(
formula=bf(
length2 | cens(censored) ~ Lmax - (Lmax - Lmin) * exp(-(k*hour)^delta), #Weibull model
Lmax + Lmin ~ (1|day),
delta + k ~ trt + (1|id), #
nl=TRUE ),
data=style,
prior=prior6,
sample_prior = "yes",
inits=inits6,
cores=4,
chains=4,
iter=2000)
Warning: It appears as if you have specified a lower bounded prior on a parameter that has no natural lower bound.
If this is really what you want, please specify argument 'lb' of 'set_prior' appropriately.
Warning occurred for prior
b_k ~ exponential(1)
b_delta ~ exponential(0.5)
Compiling Stan program...
Start sampling
starting worker pid=71112 on localhost:11685 at 09:46:45.970
starting worker pid=71127 on localhost:11685 at 09:46:46.356
starting worker pid=71141 on localhost:11685 at 09:46:46.651
starting worker pid=71155 on localhost:11685 at 09:46:46.913
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.00154 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 15.4 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 0.001593 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 15.93 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.001645 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 16.45 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4: Rejecting initial value:
Chain 4: Log probability evaluates to log(0), i.e. negative infinity.
Chain 4: Stan can't start sampling from this initial value.
Chain 4:
Chain 4: Gradient evaluation took 0.003294 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 32.94 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 215.768 seconds (Warm-up)
Chain 1: 135.726 seconds (Sampling)
Chain 1: 351.494 seconds (Total)
Chain 1:
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 267.888 seconds (Warm-up)
Chain 2: 97.123 seconds (Sampling)
Chain 2: 365.011 seconds (Total)
Chain 2:
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 270.062 seconds (Warm-up)
Chain 3: 102.446 seconds (Sampling)
Chain 3: 372.508 seconds (Total)
Chain 3:
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 265.939 seconds (Warm-up)
Chain 4: 125.089 seconds (Sampling)
Chain 4: 391.028 seconds (Total)
Chain 4:
)
Error: unexpected ')' in ")"
summary(fit_6c)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: length2 | cens(censored) ~ Lmax - (Lmax - Lmin) * exp(-(k * hour)^delta)
Lmax ~ (1 | day)
Lmin ~ (1 | day)
delta ~ trt + (1 | id)
k ~ trt + (1 | id)
Data: style (Number of observations: 1362)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Group-Level Effects:
~day (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Lmax_Intercept) 1.73 0.46 1.08 2.83 1.00 784 1330
sd(Lmin_Intercept) 0.51 0.16 0.30 0.90 1.00 1183 1699
~id (Number of levels: 72)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(delta_Intercept) 1.22 0.22 0.83 1.69 1.00 925 1704
sd(k_Intercept) 0.03 0.00 0.02 0.03 1.00 1146 2011
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Lmax_Intercept 11.44 0.23 10.97 11.89 1.00 1255 1750
Lmin_Intercept 7.47 0.15 7.19 7.77 1.00 853 1393
delta_Intercept 4.53 0.28 4.01 5.12 1.00 1021 1731
delta_trtWest -0.08 0.21 -0.50 0.33 1.00 2637 2562
delta_trtWestHeat 0.10 0.21 -0.31 0.51 1.00 2314 3031
k_Intercept 0.34 0.01 0.33 0.35 1.00 849 1459
k_trtWest -0.10 0.01 -0.12 -0.08 1.00 763 1477
k_trtWestHeat -0.03 0.01 -0.05 -0.01 1.00 785 1277
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.39 0.01 0.38 0.41 1.00 3236 3077
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 6498485 347.1 9991201 533.6 NA 9991201 533.6
Vcells 37740970 288.0 318354307 2428.9 16384 428446769 3268.8
plot(fit_6c, ask = FALSE)
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 6703482 358.1 9991201 533.6 NA 9991201 533.6
Vcells 39084255 298.2 203746757 1554.5 16384 428446769 3268.8
plot(conditional_effects(fit_6c), ask=FALSE)
pp_check(fit_6c)
Using 10 posterior samples for ppc type 'dens_overlay' by default.
Warning: Censored responses are not shown in 'pp_check'.
pred_6c <- predict(fit_6c)
pred_6c %>%
cbind(style) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_smooth(span=1) +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_smooth(aes(y=Estimate), lty=2, span=1)
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
pred_6c %>%
cbind(style) %>%
ggplot(aes(x=hour, y=length, color=trt)) +
geom_smooth(span=1) +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_smooth(aes(y=Estimate), lty=2, span=1) +
facet_wrap(~day, scales = "free_y")
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
Warning: Removed 39 rows containing non-finite values (stat_smooth).
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
ggsave("../output/style_fit_plot.pdf", width=6.5, height = 4)
Warning: Removed 13 rows containing missing values (geom_point).
save.image("../output/StyleModelFits.Rdata")
loo_compare(fit1c, fit2c, fit3c, fit_4c, fit_5c, fit_6c)
elpd_diff se_diff
fit_6c 0.0 0.0
fit_5c -157.6 20.1
fit_4c -266.4 22.7
fit3c -289.1 23.6
fit2c -443.0 33.6
fit1c -1255.8 43.9
Get the posterior
post_6c <- posterior_samples(fit_6c, pars="b.*")
posterior for inflection point
# inflection point equation
infl.pt <- function(delta, k) {
(1/k) *
((delta -1) / delta)^(1/delta)
}
post_6c$infl_intercept <- infl.pt(post_6c$b_delta_Intercept, post_6c$b_k_Intercept)
post_6c$infl_trtWest <-
infl.pt(
post_6c$b_delta_Intercept + post_6c$b_delta_trtWest,
post_6c$b_k_Intercept + post_6c$b_k_trtWest)
post_6c$infl_trtWestHeat <-
infl.pt(
post_6c$b_delta_Intercept+post_6c$b_delta_trtWestHeat,
post_6c$b_k_Intercept + post_6c$b_k_trtWestHeat)
apply(post_6c, 2, function(x) mean(x))
b_Lmax_Intercept b_Lmin_Intercept b_delta_Intercept b_delta_trtWest
1.143693e+01 7.474550e+00 4.534387e+00 -8.308510e-02
b_delta_trtWestHeat b_k_Intercept b_k_trtWest b_k_trtWestHeat
9.522676e-02 3.400632e-01 -1.007582e-01 -3.123653e-02
prior_b_Lmax prior_b_Lmin prior_b_delta_Intercept prior_b_delta_trtWest
1.149030e+01 7.500043e+00 1.967551e+00 -2.220422e-03
prior_b_delta_trtWestHeat prior_b_k_Intercept prior_b_k_trtWest prior_b_k_trtWestHeat
1.410824e-03 1.022834e+00 -6.267506e-05 -5.693913e-04
infl_intercept infl_trtWest infl_trtWestHeat
2.782406e+00 3.945179e+00 3.071197e+00
posterior probability that inflection point for trtWestHeat is greater than inflection point for East
sum( (post_6c$infl_trtWestHeat > post_6c$infl_intercept )) / nrow(post_6c)
[1] 0.99975
posterior probability that inflection point for trtWest is greater than inflection point for East
sum( (post_6c$infl_trtWest > post_6c$infl_intercept )) / nrow(post_6c)
[1] 1
posterior probability that inflection point for trtWest is greater than inflection point for trtWestHeat
sum((post_6c$infl_trtWest > post_6c$infl_trtWestHeat)) / nrow(post_6c)
[1] 1
posterior probability that k for trtWestHeat is less than East
hypothesis(fit_6c, "k_trtWestHeat < 0")
Hypothesis Tests for class b:
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.
posterior probability that k for trtWestHeat is less than East
sum( (post_6c$b_k_trtWest < 0 )) / nrow(post_6c)
[1] 1
sum( (post_6c$b_k_trtWest + post_6c$b_k_Intercept < post_6c$b_k_Intercept )) / nrow(post_6c)
[1] 1
k
post_6c %>%
as.data.frame() %>%
select(starts_with("b_k")) %>%
rename_with(~str_remove(.,"b_")) %>%
rename_with(~str_remove(.,"trt")) %>%
rename(k_East=k_Intercept) %>%
mutate(k_West=k_West+k_East,
k_WestHeat=k_WestHeat+k_East) %>%
summarize(across(.fns = list(mean = mean,
lower = ~ rethinking::HPDI(., 0.95)[1],
upper = ~ rethinking::HPDI(., 0.95)[2])
)) %>%
pivot_longer(cols=everything(), names_to=c("trt", "parameter"), names_pattern = "k_(.*)_(.*)") %>%
pivot_wider(names_from = "parameter", values_from = "value") %>%
ggplot(aes(x=trt, y=mean, ymin=lower, ymax=upper)) +
geom_col(fill="skyblue") +
geom_errorbar(width=.5) +
ylab("k") +
xlab("Condition") +
theme_bw()
ggsave("../output/style_k_plot.pdf", height=3, width = 3)
inflection
post_6c %>%
as.data.frame() %>%
select(starts_with("infl")) %>%
rename_with(~str_remove(.,"trt")) %>%
rename_with(~str_remove(.,"infl_")) %>%
rename(East=intercept) %>%
summarize(across(.fns = list(mean = mean,
lower = ~ rethinking::HPDI(., 0.95)[1],
upper = ~ rethinking::HPDI(., 0.95)[2])
)) %>%
pivot_longer(cols=everything(), names_to=c("trt", "parameter"), names_pattern = "(.*)_(.*)") %>%
pivot_wider(names_from = "parameter", values_from = "value") %>%
ggplot(aes(x=trt, y=mean, ymin=lower, ymax=upper)) +
geom_col(fill="red") +
geom_errorbar(width=.5) +
ylab("Inflection Point") +
xlab("Condition") +
theme_bw()
ggsave("../output/style_inflection_plot.pdf", height=3, width = 3)
posterior probability that k for trtWest is less than k for trtWestHEat
sum( (post_6c$b_k_trtWest < post_6c$b_k_trtWestHeat )) / nrow(post_6c)
[1] 1